### My script starts here.  This script will track temperature anomalies through time across all cells in an ASCII
### The user defined threshold will be 28 degrees.  This script will calculate for all cells the longest streak above threshold
### And total days above threshold

# First define directories and install any packages that may be necessary

in.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/SpatialTMax/daysabove/'
out.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/SpatialTMax/daysabove/'
setwd(in.dir)

#The following function was written by JJV
#It will identify numbers in a vector list above a user defined threshhold and group them into unique clusters

extreme.clust = function(x,y,thresh.type){
if (thresh.type=='>=') z = which(x>=y)
if (thresh.type=='>') z = which(x>y)
if (thresh.type=='<=') z = which(x<=y)
if (thresh.type=='<') z = which(x<y)
if (thresh.type=='==') z = which(x==y)
if (length(z)>1){
z.diff = diff(z)
z.clust = NULL; tclust=1
for (i in 1:length(z.diff)){
z.clust = c(z.clust,tclust)
if (z.diff[i]>1) tclust = tclust + 1
if (i == length(z.diff)) z.clust = c(z.clust,tclust)
}
out = rep(NA,length(x))
out[z] = z.clust
return(out)
} else if (length(z)==1) {
out = rep(NA,length(x))
out[z] = 1
return(out)
} else {
return(NULL)
}
}

clust.data = read.csv('clusterdata5.csv', header=T) # Read in map data in expanded format
threshold = 280 #Set threshold (280 not 28 because all values were multiplied by 10 and rounded to become integer)
pos = 1:365  #Number of days to track pattern across i.e. 1:365

# Create a matrix that has the same amount of rows as clust.data and 4 columns (row/col/streak/total

t.mat = matrix(0,nr=nrow(clust.data),nc=4)

# Populate matrix with row column values from clust.data

t.mat[,1]=clust.data[,1]
t.mat[,2]=clust.data[,2]

for (i in 1:nrow(clust.data))

    {

    cells.above = extreme.clust(clust.data[i,5:369],threshold,'>=')

    cat(i,'...\n')

    if(is.null(cells.above))

      {

      streak.above = 0

      total.above = 0

      t.mat[i,3]= streak.above[1]

      t.mat[i,4]= total.above[1]

      }

        else

        {

        out=aggregate(pos,list(cluster=cells.above),length)

        streak.above = max(out[,2])

        total.above = sum(out[,2])

        t.mat[i,3]= streak.above[1]

        t.mat[i,4]= total.above[1]

        }

  }

write.csv(x=t.mat, file=paste(out.dir,'summary5.csv',sep=""),row.names=F)